/*
9.27
1.function matrix_calc_taskA,matrix_calc_taskB
2."N" means iteration times
10.2
1.merge the for loop in taskB
2.move memory access variable out of "for loop".

10.19
4 row at the same time 
*/
#ifndef _CALC_H_
#define _CALC_H_
#include<immintrin.h>
#ifdef __cplusplus
extern "C"{	
#endif

typedef struct {
	int* 			rowArray;
	const int* 		rowOffset;
	int 			rowArraySize;
	const int* 		columnIndice;
	const double* 	S;
	const double*	valueNormalMatrix;
	double* 		Id;
}TaskMatrixInfoA;

typedef struct {
	const double * 	valueSpiceMatrix;//N*2
	const int* 		rowOffset; 		//n+1
	const int* 		columnIndice; 	//N

	double* 		A;	//N
	double*		 	S;	//nx1
	double*			R;	//nx1
	double* 		H;	//nx1
	const double* 	D;	//N*2
	double* 		IC;	//nx1
	double* 		IG;	//nx1
	double 		alpha;
	int * 			rowArray; //kx1
	int 			rowArraySize;
	void* 			hdl;  //? what's this?
}TaskMatrixInfoB;

void matrix_calc_taskA(TaskMatrixInfoA** listDataList,int N)
{
	for (int i = 0 ;i < N; ++i)
	{	
		int row_array_size= listDataList[i]->rowArraySize;	
		TaskMatrixInfoA * list_Data_List=listDataList[i];
		const double *vNM = list_Data_List->valueNormalMatrix;
		const double *S = list_Data_List->S;
		const int *cI = list_Data_List->columnIndice;
		const int *rOf= list_Data_List->rowOffset;
		int *rA =list_Data_List->rowArray;
		double *Id = list_Data_List->Id;

//		int left = row_array_size % 8;
		int left = row_array_size &0x07;
		int it;

		for ( it = 0;it<left;it++){
			const int node = list_Data_List->rowArray[it];
			int colum_num = list_Data_List->rowOffset[node + 1];
			double &id=list_Data_List->Id[node];

			for (int j = list_Data_List->rowOffset[node];j < colum_num;++j){   
				id += vNM[j]*S[cI[j]];
			}
		}

		double tmp[8] __attribute__((aligned(64)));
		__m512d m1,m2,m3;
		int node8[8];
		int jstart[8];
		int jend[8];

		for (it = left ; it < row_array_size; it=it + 8)
		{
			node8[0]=rA[it];  node8[1]=rA[it+1];
			node8[2]=rA[it+2];node8[3]=rA[it+3];
			node8[4]=rA[it+4];node8[5]=rA[it+5];
			node8[6]=rA[it+6];node8[7]=rA[it+7];
			jstart[0]=rOf[node8[0]];jend[0]=rOf[node8[0]+1];
			jstart[1]=rOf[node8[1]];jend[1]=rOf[node8[1]+1];
			jstart[2]=rOf[node8[2]];jend[2]=rOf[node8[2]+1];
			jstart[3]=rOf[node8[3]];jend[3]=rOf[node8[3]+1];
			jstart[4]=rOf[node8[4]];jend[4]=rOf[node8[4]+1];
			jstart[5]=rOf[node8[5]];jend[5]=rOf[node8[5]+1];
			jstart[6]=rOf[node8[6]];jend[6]=rOf[node8[6]+1];
			jstart[7]=rOf[node8[7]];jend[7]=rOf[node8[7]+1];

			m1 = _mm512_set_pd(vNM[jstart[7]],vNM[jstart[6]],vNM[jstart[5]],vNM[jstart[4]],vNM[jstart[3]],vNM[jstart[2]],vNM[jstart[1]],vNM[jstart[0]]);
			m2 = _mm512_set_pd(S[cI[jstart[7]]],S[cI[jstart[6]]],S[cI[jstart[5]]],S[cI[jstart[4]]],S[cI[jstart[3]]],S[cI[jstart[2]]],S[cI[jstart[1]]],S[cI[jstart[0]]]);
			m3 = _mm512_mul_pd(m1,m2);
			_mm512_store_pd(tmp,m3);

			for (int it2=0;it2<8;it2++){
				for(int tmpj = jstart[it2]+1;tmpj<jend[it2];tmpj++){
					tmp[it2]+=vNM[tmpj]*S[cI[tmpj]];	
				}
			}
			Id[node8[0]]+=tmp[0];	Id[node8[1]]+=tmp[1];
			Id[node8[2]]+=tmp[2];	Id[node8[3]]+=tmp[3];
			Id[node8[4]]+=tmp[4];	Id[node8[5]]+=tmp[5];
			Id[node8[6]]+=tmp[6];	Id[node8[7]]+=tmp[7];
		}
	}
}

//Task B interface
void matrix_calc_taskB(TaskMatrixInfoB** listDataList,int N)
{
	for (int i = 0;i<N;i++){
// 定义临时变量
		TaskMatrixInfoB * list_Data_List = listDataList[i];
		int row_array_size = list_Data_List->rowArraySize;
		const double alpha = list_Data_List->alpha;
		double* A = list_Data_List->A;
		const double* S = list_Data_List->S;
		const int* cI=list_Data_List->columnIndice;
		const double* vSM=list_Data_List->valueSpiceMatrix;
		int *rA = list_Data_List->rowArray;
		const int* rOf = list_Data_List->rowOffset;

// 第二层循环
		int left = row_array_size & 0x07;
		int it;

//计算前left行，剩下的行数为8的倍数
		for (it =0;it<left;it++){
			int node = rA[it];
			const int k1 = node << 1;
			int jend = rOf[node+1];
			double ig_tmp=0;
			double ic_tmp=0;

			for (int p = rOf[node];p < jend;++p){
				const int k = p << 1;
				double cond =vSM[k];
				double cap  =vSM[k+1];					
				ig_tmp += cond * S[cI[p]];
				ic_tmp += cap * S[cI[p]];
				//(9)
				A[p] = cond +  cap * alpha;
			}
			list_Data_List->IG[node] +=ig_tmp; 
 			list_Data_List->IC[node] +=ic_tmp; 
			list_Data_List->R [node] = list_Data_List->D[k1] -ig_tmp;
			list_Data_List->H [node] = list_Data_List->D[k1+1]-ic_tmp;
		}
//计算后面的8的倍数的行数		
		double ig_tmp[8] __attribute__((aligned(64)));
		double ic_tmp[8] __attribute__((aligned(64)));
		__m512d mg,mc,ms,m3;
		int node8[8];
		int jstart[8];
		int jend[8];
		int k1[8];

		for (it = left;it<row_array_size;it+=8){
			node8[0]=rA[it];  node8[1]=rA[it+1];
			node8[2]=rA[it+2];node8[3]=rA[it+3];
			node8[4]=rA[it+4];node8[5]=rA[it+5];
			node8[6]=rA[it+6];node8[7]=rA[it+7];
			jstart[0]=rOf[node8[0]];jend[0]=rOf[node8[0]+1];
			jstart[1]=rOf[node8[1]];jend[1]=rOf[node8[1]+1];
			jstart[2]=rOf[node8[2]];jend[2]=rOf[node8[2]+1];
			jstart[3]=rOf[node8[3]];jend[3]=rOf[node8[3]+1];
			jstart[4]=rOf[node8[4]];jend[4]=rOf[node8[4]+1];
			jstart[5]=rOf[node8[5]];jend[5]=rOf[node8[5]+1];
			jstart[6]=rOf[node8[6]];jend[6]=rOf[node8[6]+1];
			jstart[7]=rOf[node8[7]];jend[7]=rOf[node8[7]+1];
			k1[0]= node8[0]<<1;k1[1]= node8[1]<<1;
			k1[2]= node8[2]<<1;k1[3]= node8[3]<<1;
			k1[4]= node8[4]<<1;k1[5]= node8[5]<<1;
			k1[6]= node8[6]<<1;k1[7]= node8[7]<<1;
//计算每一行的第一个乘积
			mg = _mm512_set_pd(
				vSM[jstart[7]<<1],vSM[jstart[6]<<1],
				vSM[jstart[5]<<1],vSM[jstart[4]<<1],
				vSM[jstart[3]<<1],vSM[jstart[2]<<1],
				vSM[jstart[1]<<1],vSM[jstart[0]<<1]);
			mc = _mm512_set_pd(
				vSM[(jstart[7]<<1)+1],vSM[(jstart[6]<<1)+1],
				vSM[(jstart[5]<<1)+1],vSM[(jstart[4]<<1)+1],
				vSM[(jstart[3]<<1)+1],vSM[(jstart[2]<<1)+1],
				vSM[(jstart[1]<<1)+1],vSM[(jstart[0]<<1)+1]);
			ms = _mm512_set_pd(
				S[cI[jstart[7]]],S[cI[jstart[6]]],
				S[cI[jstart[5]]],S[cI[jstart[4]]],
				S[cI[jstart[3]]],S[cI[jstart[2]]],
				S[cI[jstart[1]]],S[cI[jstart[0]]]);
			m3 = _mm512_mul_pd(mg,ms);
			_mm512_store_pd(ig_tmp,m3);
			m3 = _mm512_mul_pd(mc,ms);
			_mm512_store_pd(ic_tmp,m3);
//计算每行剩下的乘积
			A[jstart[0]]=vSM[jstart[0]<<1]+vSM[(jstart[0]<<1)+1]*alpha;
			A[jstart[1]]=vSM[jstart[1]<<1]+vSM[(jstart[1]<<1)+1]*alpha;
			A[jstart[2]]=vSM[jstart[2]<<1]+vSM[(jstart[2]<<1)+1]*alpha;
			A[jstart[3]]=vSM[jstart[3]<<1]+vSM[(jstart[3]<<1)+1]*alpha;
			A[jstart[4]]=vSM[jstart[4]<<1]+vSM[(jstart[4]<<1)+1]*alpha;
			A[jstart[5]]=vSM[jstart[5]<<1]+vSM[(jstart[5]<<1)+1]*alpha;
			A[jstart[6]]=vSM[jstart[6]<<1]+vSM[(jstart[6]<<1)+1]*alpha;
			A[jstart[7]]=vSM[jstart[7]<<1]+vSM[(jstart[7]<<1)+1]*alpha;

			for(int it2=0;it2<8;it2++){
				for(int tmpj = jstart[it2]+1;tmpj<jend[it2];tmpj++){
					ig_tmp[it2]+= vSM[tmpj<<1]*S[cI[tmpj]];
					ic_tmp[it2]+= vSM[(tmpj<<1)+1]*S[cI[tmpj]];
					A[tmpj]  =vSM[tmpj<<1]+vSM[(tmpj<<1)+1]*alpha; 
				}
			}
			list_Data_List->IG[node8[0]]+=ig_tmp[0];
			list_Data_List->IG[node8[1]]+=ig_tmp[1];
			list_Data_List->IG[node8[2]]+=ig_tmp[2];
			list_Data_List->IG[node8[3]]+=ig_tmp[3];
			list_Data_List->IG[node8[4]]+=ig_tmp[4];
			list_Data_List->IG[node8[5]]+=ig_tmp[5];
			list_Data_List->IG[node8[6]]+=ig_tmp[6];
			list_Data_List->IG[node8[7]]+=ig_tmp[7];

			list_Data_List->IC[node8[0]]+=ic_tmp[0];
			list_Data_List->IC[node8[1]]+=ic_tmp[1];
			list_Data_List->IC[node8[2]]+=ic_tmp[2];
			list_Data_List->IC[node8[3]]+=ic_tmp[3];
			list_Data_List->IC[node8[4]]+=ic_tmp[4];
			list_Data_List->IC[node8[5]]+=ic_tmp[5];
			list_Data_List->IC[node8[6]]+=ic_tmp[6];
			list_Data_List->IC[node8[7]]+=ic_tmp[7];	
			
			list_Data_List->R[node8[0]] =list_Data_List->D[k1[0]]-ig_tmp[0];
			list_Data_List->R[node8[1]] =list_Data_List->D[k1[1]]-ig_tmp[1];
			list_Data_List->R[node8[2]] =list_Data_List->D[k1[2]]-ig_tmp[2];
			list_Data_List->R[node8[3]] =list_Data_List->D[k1[3]]-ig_tmp[3];
			list_Data_List->R[node8[4]] =list_Data_List->D[k1[4]]-ig_tmp[4];
			list_Data_List->R[node8[5]] =list_Data_List->D[k1[5]]-ig_tmp[5];
			list_Data_List->R[node8[6]] =list_Data_List->D[k1[6]]-ig_tmp[6];
			list_Data_List->R[node8[7]] =list_Data_List->D[k1[7]]-ig_tmp[7];

			list_Data_List->H[node8[0]] =list_Data_List->D[k1[0]+1]-ic_tmp[0];
			list_Data_List->H[node8[1]] =list_Data_List->D[k1[1]+1]-ic_tmp[1];
			list_Data_List->H[node8[2]] =list_Data_List->D[k1[2]+1]-ic_tmp[2];
			list_Data_List->H[node8[3]] =list_Data_List->D[k1[3]+1]-ic_tmp[3];
			list_Data_List->H[node8[4]] =list_Data_List->D[k1[4]+1]-ic_tmp[4];
			list_Data_List->H[node8[5]] =list_Data_List->D[k1[5]+1]-ic_tmp[5];
			list_Data_List->H[node8[6]] =list_Data_List->D[k1[6]+1]-ic_tmp[6];
			list_Data_List->H[node8[7]] =list_Data_List->D[k1[7]+1]-ic_tmp[7];			
		}
/*
		for(int it = 0;it < row_array_size;++it)
		{
			int row = list_Data_List->rowArray[it];
			const int k1 = row * 2;
			int colum_num = list_Data_List->rowOffset[row + 1];
							
			double ig_tmp=0;
			double ic_tmp=0;
			for (int p = list_Data_List->rowOffset[row];p < colum_num;++p){
				const int k = p << 1;
				double cond =valueSpiceMatrix[k];
				double cap  =valueSpiceMatrix[k+1];					
				ig_tmp += cond * S[columnIndice[p]];
				ic_tmp += cap * S[columnIndice[p]];
				//(9)
				A[p] = cond +  cap * alpha;
			}
 			//(7)(8)
 			list_Data_List->IG[row] +=ig_tmp; 
 			list_Data_List->IC[row] +=ic_tmp; 
			list_Data_List->R[row] = list_Data_List->D[k1]-ig_tmp;
			list_Data_List->H[row] = list_Data_List->D[k1+1]-ic_tmp;
		} 
*/
	}
}

#ifdef __cplusplus
}
#endif

#endif
